Financial Econometrics I: Project - Pavlína Křenková and Bernhard Brunner¶
June 2024
# Packages
shhh = suppressPackageStartupMessages #load quietly
shhh(library(zip))
shhh(library(ggplot2))
shhh(library(gridExtra))
shhh(library(forecast))
shhh(library(xts))
shhh(library(tseries))
shhh(library(highfrequency))
shhh(library(lmtest))
shhh(library(sandwich))
shhh(library(zoo))
shhh(library(rugarch))
shhh(library(DescTools))
shhh(library(car))
1) Introduction¶
This project estimates and compares several econometric models analyzing the volatility of the American Tower Corp stock returns. 3 main sections are included: data description, focusing on exploration of the data using summary statistics and simple plots, in-sample fit section in which 6 different models are fitted and evaluated, and lastly the forecasts section where we compare the out-of-sample forecasting performance of all the models using both the rolling and expanding window forecasting schemes.
2) Data description¶
Load the data into R¶
unzip("data_project_2024.zip", files = "30.RData", exdir = "temp_data")
# Load the .RData file
load("temp_data/30.RData")
# Return all variables and functions defined in current working directory
ls()
- 'amt'
- 'shhh'
Descriptive statistics¶
head(amt)
summary(amt)
nrow(amt) # Get number of observations
typeof(amt) # Get data type of time series object
ret RV RV_p RV_n RS
2010-01-05 0.0002304634 0.011009479 0.009259483 0.005955720 1.350981427
2010-01-06 0.0168678327 0.011202805 0.008230314 0.007600314 -0.216028901
2010-01-07 0.0042856864 0.011699048 0.008224708 0.008319971 -0.006426756
2010-01-08 -0.0011261662 0.009363287 0.007205883 0.005978829 0.753857266
2010-01-11 0.0109791088 0.011245783 0.006254063 0.009346354 -0.759894442
2010-01-12 -0.0252775180 0.011362306 0.006760609 0.009132151 -0.795626739
RK
2010-01-05 5.926201
2010-01-06 6.845365
2010-01-07 4.509419
2010-01-08 6.858052
2010-01-11 5.162329
2010-01-12 4.448177
Index ret RV
Min. :2010-01-05 Min. :-0.1035571 Min. :0.003823
1st Qu.:2011-07-04 1st Qu.:-0.0071153 1st Qu.:0.007748
Median :2013-01-10 Median : 0.0006953 Median :0.009442
Mean :2013-01-08 Mean : 0.0005369 Mean :0.010519
3rd Qu.:2014-07-14 3rd Qu.: 0.0079566 3rd Qu.:0.011820
Max. :2016-01-22 Max. : 0.0701801 Max. :0.055211
RV_p RV_n RS RK
Min. :0.002432 Min. :0.002027 Min. :-5.71549 Min. : 2.054
1st Qu.:0.005315 1st Qu.:0.005263 1st Qu.:-0.59881 1st Qu.: 3.720
Median :0.006521 Median :0.006564 Median : 0.01357 Median : 4.956
Mean :0.007430 Mean :0.007262 Mean : 0.04922 Mean : 6.633
3rd Qu.:0.008385 3rd Qu.:0.008277 3rd Qu.: 0.67202 3rd Qu.: 7.302
Max. :0.054001 Max. :0.039343 Max. : 7.61326 Max. :64.104
By calling the head() function we display the first five rows of the American Tower Cooperation (AMT) time series object . The displayed variables are Returns, Realized Volatility, Positive Realized Semi-Volatility, Negative Realized Semi-Volatility, Realized Skewness, Realized Kurtosis.
By calling the summary() function we obtain the minimum, maximumum, mean, median, first-, and third quantile of each variable in the times series object.
Our data spans from Jan 2010 to Jan 2016. Overall, the values seem reasonable, with the mean of the returns being slightly positive, but around zero.
Plot the data¶
plot_dta = function(column, title) {
ggplot(amt, aes(x = index(amt), y = get(column))) +
geom_line(color = "#194f80") +
labs(title = title, x = "Time", y = column) +
theme_minimal()
}
# Assign variables and column names
plots = list()
titles = c("Returns", "Realized Volatility", "Positive Realized Semi-Volatility",
"Negative Realized Semi-Volatility", "Realized Skewness", "Realized Kurtosis")
# Loop over each variables to create time series plots iteratively
for (i in seq_along(names(amt))) {
plots[[i]] = plot_dta(names(amt)[i], titles[i])
}
options(repr.plot.width = 17, repr.plot.height = 10)
grid.arrange(grobs = plots, ncol = 2, nrow = 3)
We plotted the Returns, Realized Volatility, Positive Realized Semi-Volatility, Negative Realized Semi-Volatility, Realized Skewness, and Realized Kurtosis above. Plotting the Returns indicates a mean around zero and changes in volatility over time with high-volatility clusters around mid-2010, mid-2011, mid-2013, and mid-2015. These high-volatility clusters are also shown by the Realized Volatility plot, where the Realized Volatility peaks roughly coincide with the high-volatility clusters in the Returns plot, the RV plot also shows several jumps. Positive and Negative Semi-Realized Volatility seem to be closely related to the Realized Volatility plot, due to their mathematical definition of either only considering positive or negative returns. Realized Skewness and Kurtosis also vary over time.
options(repr.plot.width = 10, repr.plot.height = 5)
ggplot(amt, aes(x = ret)) +
geom_histogram(aes(y = after_stat(density)), fill = "#194f80", bins = 30, color = "black") +
geom_density(color = "red", linewidth = 1) +
stat_function(fun = dnorm, args = list(mean = mean(amt$ret, na.rm = TRUE), sd = sd(amt$ret, na.rm = TRUE)), color = "green", linewidth = 1) +
ggtitle("Histogram of amt_ret with Normal Distribution Curve") +
theme_minimal()
The green curve is a normal distribution pdf that would correspond to the sample mean and sd of the of the returns. Clearly, the red line, which estimates the return density, exerts a much higher kurtosis, commonly observed with stock returns.
Plot the autocorrelation (ACF) and partial-autocorrelation (PACF) functions¶
# Plot ACF and PACF functions for Returns
options(repr.plot.width = 10, repr.plot.height = 5)
par(mfrow = c(1, 2))
Acf(amt$ret, main="ACF for Returns")
Pacf(amt$ret, main="PACF for Returns")
adf.test(amt$ret)
Warning message in adf.test(amt$ret): “p-value smaller than printed p-value”
Augmented Dickey-Fuller Test data: amt$ret Dickey-Fuller = -12.231, Lag order = 11, p-value = 0.01 alternative hypothesis: stationary
The ACF show four statistically significant autocorrelations at lag 1, 3, 18, and 27. The PACF shows seven statistically significant partial autocorrelations at lag 1, 2, 3, 5, 18, 27, and 30. This suggests that the return series is not white noise and there are some significant (partial) autocorrelations, indicating that the series has some predictability based on past values. The $H_0$ of the adf test is rejected, suggesting the series is stationary.
# Plot ACF, and PACF of Realized Volatility variables
par(mfrow = c(3, 2))
options(repr.plot.width = 10, repr.plot.height = 10)
Acf(amt$RV, main="ACF for Realized Volatility")
Pacf(amt$RV, main="PACF for Realized Volatility")
Acf(amt$RV_p, main="ACF for Positive Realized semi-Volatility")
Pacf(amt$RV_p, main="PACF for Positive Realized semi-Volatility")
Acf(amt$RV_n, main="ACF for Negative Realized semi-Volatility")
Pacf(amt$RV_n, main="PACF for Negative Realized semi-Volatility")
Realized Volatility, Positive Realized Semi-Volatility, and Negative Realized Semi-Volatility have similar ACF and PACF functions. The ACF for all three variables is slowly decaying, which could possibly indicate non-stationarity; however we do reject the H_0 of the adf.test below. Furthermore, the PACFs show statistically significant partial autocorrelations, but the decay is much faster, indicating that the series could possibly follow some AR process.
adf.test(amt$RV)
adf.test(amt$RV_n)
adf.test(amt$RV_p)
Warning message in adf.test(amt$RV): “p-value smaller than printed p-value”
Augmented Dickey-Fuller Test data: amt$RV Dickey-Fuller = -6.5222, Lag order = 11, p-value = 0.01 alternative hypothesis: stationary
Warning message in adf.test(amt$RV_n): “p-value smaller than printed p-value”
Augmented Dickey-Fuller Test data: amt$RV_n Dickey-Fuller = -6.7825, Lag order = 11, p-value = 0.01 alternative hypothesis: stationary
Warning message in adf.test(amt$RV_p): “p-value smaller than printed p-value”
Augmented Dickey-Fuller Test data: amt$RV_p Dickey-Fuller = -6.8486, Lag order = 11, p-value = 0.01 alternative hypothesis: stationary
3) In-sample fit¶
(a) AR(1)-RV¶
where:
- $RV_t$ is the realized volatility at time ( t ),
- $\phi_0$ is the intercept,
- $ \phi_1 $ is the coefficient of lag 1 of realized volatility,
- $ \epsilon_t $ is the error term, assumed to be normally distributed white noise.
# Fit the AR(1)-RV model
rv_ar1_fit = Arima(amt$RV, order=c(1,0,0))
rv_ar1_fit
Series: amt$RV
ARIMA(1,0,0) with non-zero mean
Coefficients:
ar1 mean
0.456 0.0105
s.e. 0.023 0.0002
sigma^2 = 1.759e-05: log likelihood = 6083.53
AIC=-12161.07 AICc=-12161.05 BIC=-12145.13
Our estimated coefficient for lag 1 Realized Volatility is 0.456 with a standard error of 0.023, indicating a positive relationship between the lagged and current values. The estimated intercept is 0.015 with a standard error of 0.0002. Both estimates are statistically significant.
# Plot the true versus AR(1)-RV fitted values
options(repr.plot.width = 10, repr.plot.height = 5)
amt$rv_ar1_fit_val = xts(fitted(rv_ar1_fit), order.by = index(amt))
head(amt$rv_ar1_fit_val)
plot_rv_ar = plot(amt$RV, type = 'l', col = "#194f80", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - AR(1)-RV")
plot_rv_ar = lines(amt$rv_ar1_fit_val, col = "green", lwd = 1)
plot_rv_ar = addLegend("topright", legend.names = c("True Values", "Fitted Values"), col = c("#194f80", "green"), lty = 1, lwd = 2)
plot_rv_ar
rv_ar1_fit_val 2010-01-05 0.010573280 2010-01-06 0.010742866 2010-01-07 0.010831030 2010-01-08 0.011057337 2010-01-11 0.009992136 2010-01-12 0.010850630
Here, we plot the true values against the AR(1)-RV fitted values. The graph indicates a discrepancy between the true and fitted values. Specifically, it seems that the AR(1)-RV model overestimates the volatility, while not capturing the jumps and the more detailed structure of the realized volatility; however it is able to capture the general pattern and approximate magnitude.
(b) Heterogeneous Autoregressive Model¶
where:
- $RV_t$ is the daily realized volatility at time ( t ),
- $RV_{t}^{(5)}$ is the weekly realized volatility at time ( t ),
- $RV_{t}^{(22)}$ is the monthly realized volatility at time ( t ),
- $\alpha$ is the intercept,
- $ \beta_1 $ is the coefficient of lag 1 of daily realized volatility,
- $ \beta_2 $ is the coefficient of lag 1 of weekly realized volatility,
- $ \beta_3 $ is the coefficient of lag 1 of monthly realized volatility,
- $ u_t $ is the error term, assumed to be normally distributed white noise.
# Fit the HAR model and print estimates
HAR_fit = HARmodel(data = amt$RV , periods = c(1,5,22))
summary(HAR_fit)
Call:
"RV1 = beta0 + beta1 * RV1 + beta2 * RV5 + beta3 * RV22"
Residuals:
Min 1Q Median 3Q Max
-0.012629 -0.002062 -0.000618 0.001081 0.045830
Coefficients:
Estimate Std. Error t value Pr(>|t|)
beta0 0.0018621 0.0003657 5.092 4.01e-07 ***
beta1 0.2298062 0.0474828 4.840 1.44e-06 ***
beta2 0.2777390 0.0796523 3.487 0.000503 ***
beta3 0.3155588 0.0893409 3.532 0.000425 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.004013 on 1474 degrees of freedom
Multiple R-squared: 0.2775, Adjusted R-squared: 0.276
F-statistic: 188.7 on 3 and 1474 DF, p-value: < 2.2e-16
# Show the estimated coefficient and significance levels manually
amt$RV_1 = lag(amt$RV)
amt$RV_5 = lag(rollapply(amt$RV, width = 5, FUN = mean, fill = NA, align = "right"))
amt$RV_22 = lag(rollapply(amt$RV, width = 22, FUN = mean, fill = NA, align = "right"))
har_manual = lm(RV ~ RV_1 + RV_5 + RV_22, data = amt)
coeftest(har_manual, vcov = NeweyWest(har_manual, lag = 22))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.00186209 0.00036572 5.0916 4.008e-07 ***
RV_1 0.22980621 0.04748279 4.8398 1.436e-06 ***
RV_5 0.27773900 0.07965227 3.4869 0.000503 ***
RV_22 0.31555884 0.08934093 3.5321 0.000425 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
All coefficients for the HAR model are statistically significant and positive, indicating that the lagged average monthly and weekly realized volatility are also useful for the modelling.
head(xts(fitted(HAR_fit), order.by = as.Date(names(fitted(HAR_fit)))))
[,1] 2010-02-05 0.01361216 2010-02-08 0.01442904 2010-02-09 0.01360447 2010-02-10 0.01328010 2010-02-11 0.01295733 2010-02-12 0.01441048
amt$HAR_fit_val = xts(fitted(HAR_fit), order.by = as.Date(names(fitted(HAR_fit))))
plot_HAR = plot(amt$RV, type = 'l', col = "#194f80", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - HAR")
plot_HAR = lines(amt$HAR_fit_val, col = "green", lwd = 1)
plot_HAR = addLegend("topright", legend.names = c("True Values", "Fitted Values"), col = c("#194f80", "green"), lty = 1, lwd = 2)
plot_HAR
Here, we plot the true versus fitted HAR values. The fit appears quite similar to before, but it seems that the HAR estimates better fits the data than the AR(1)-RV model. Specifically, it seems from the figure that the HAR model does not overestimate mean volatility as strongly as the AR(1)-RV model.
(c) Heterogeneous Autoregressive Model with Realized Semi-volatility¶
where:
- $RV_t$ is the daily realized volatility at time ( t ),
- $RS^+_{t} $ is the daily positive realized semi-volatility at time ( t ),
- $RS^-_{t}$ is the daily negative realized semi-volatility at time ( t ),
- $RV_{t}^{(5)}$ is the weekly realized volatility at time ( t ),
- $RV_{t}^{(22)}$ is the monthly realized volatility at time ( t ),
- $\alpha$ is the intercept,
- $ \beta_1 $ is the coefficient of lag 1 of daily positive realized semi-volatility,
- $ \beta_2 $ is the coefficient of lag 1 of daily positive realized semi-volatility,
- $ \beta_3 $ is the coefficient of lag 1 of weekly realized volatility,
- $ \beta_4 $ is the coefficient of lag 1 of monthly realized volatility,
- $ u_t $ is the error term, assumed to be normally distributed white noise.
amt$RV_p_l = lag(amt$RV_p)
amt$RV_n_l = lag(amt$RV_n)
har_rs_fit = lm(RV ~ RV_p_l + RV_n_l + RV_5 + RV_22, data = amt)
coeftest(har_rs_fit, vcov = NeweyWest(har_rs_fit, lag = 22))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.00172761 0.00036179 4.7752 1.974e-06 ***
RV_p_l 0.03473701 0.03482037 0.9976 0.3186342
RV_n_l 0.37126020 0.07133673 5.2043 2.222e-07 ***
RV_5 0.24759172 0.07752767 3.1936 0.0014349 **
RV_22 0.30761444 0.08464799 3.6340 0.0002886 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
linearHypothesis(har_rs_fit, "RV_p_l = RV_n_l") #
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) | |
|---|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 1474 | 0.02365038 | NA | NA | NA | NA |
| 2 | 1473 | 0.02331648 | 1 | 0.0003339057 | 21.09423 | 4.745098e-06 |
Previously we estimated a positive coefficient for the $RV_{t-1}$, while only the coefficient on the negative realized semivolatility is statistically significant here. It seems that $\hat{\beta_1} \neq \hat{\beta_2}$, hence we were able to get new information from using semi-volatilities. We can see that all estimates except the one of $\beta_1$ are statistically significant and postive. The fact that the negative semi-volatility appears to have a larger impact than the positive one is in line with the theory that the reaction to "bad news" is stronger.
# Plot the true versus fitted values of HAR with semi-realized volatility
amt$HAR_semi_fit_val = xts(fitted(har_rs_fit), order.by = as.Date(names(fitted(har_rs_fit))))
plot_HAR_semi = plot(amt$RV, type = 'l', col = "#194f80", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - HAR w/ simi-vol")
plot_HAR_semi = lines(amt$HAR_semi_fit_val, col = "green", lwd = 1)
plot_HAR_semi = addLegend("topright", legend.names = c("True Values", "Fitted Values"), col = c("#194f80", "green"), lty = 1, lwd = 2)
plot_HAR_semi
The plot appears similar to before. We will also observe the differences when we plot the fitted values of all the models in one graph.
(d) Heterogeneous Autoregressive Model with Realized Skewness and Kurtosis¶
where:
- $RV_t$ is the daily realized volatility at time ( t ),
- $RV_{t}^{(5)}$ is the weekly realized volatility at time ( t ),
- $RV_{t}^{(22)}$ is the monthly realized volatility at time ( t ),
- $RKurt_{t} $ is the daily realized kurtosis at time ( t ),
- $RSkew_{t}$ is the daily realized skewness at time ( t ),
- $\alpha$ is the intercept,
- $ \beta_1 $ is the coefficient of lag 1 of daily realized volatility,
- $ \beta_2 $ is the coefficient of lag 1 of weekly realized volatility,
- $ \beta_3 $ is the coefficient of lag 1 of monthly realized volatility,
- $ \beta_4 $ is the coefficient of lag 1 of dailty realized kurtosis,
- $ \beta_5 $ is the coefficient of lag 1 of dailty realized skewness,
- $ u_t $ is the error term, assumed to be normally distributed white noise.
# Estimate the HAR with Realized Skewness and Kurtosis
amt$RS_l = lag(amt$RS)
amt$RK_l = lag(amt$RK)
har_rsk_rkurt_fit = lm(RV ~ RV_1 + RV_5 + RV_22 + RS_l + RK_l, data = amt)
coeftest(har_rsk_rkurt_fit, vcov = NeweyWest(har_rsk_rkurt_fit, lag = 22))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.5702e-03 3.6358e-04 7.0693 2.396e-12 ***
RV_1 3.1487e-01 5.2423e-02 6.0064 2.387e-09 ***
RV_5 2.4407e-01 7.3900e-02 3.3027 0.0009804 ***
RV_22 2.6958e-01 8.6922e-02 3.1014 0.0019626 **
RS_l -1.2198e-04 8.4773e-05 -1.4389 0.1503777
RK_l -1.1444e-04 2.7186e-05 -4.2096 2.713e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
In our case, we do not have a statistically significant estimate for realized skewness. The results indicate that larger kurtosis in previous period should on average lead to lower values of $RV$. This might be a bit of an unexpected result since fatter tails mean higher risk which we would associate with higher volatility.
# Plot the true versus fitted values of the HAR with Realized Skewness and Kurtosis
amt$HAR_rs_rk_fit_val = xts(fitted(har_rsk_rkurt_fit), order.by = as.Date(names(fitted(har_rsk_rkurt_fit))))
plot_HAR_rk_rs = plot(amt$RV, type = 'l', col = "#194f80", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - HAR w/ RK & RS")
plot_HAR_rk_rs = lines(amt$HAR_rs_rk_fit_val, col = "green", lwd = 1)
plot_HAR_rk_rs = addLegend("topright", legend.names = c("True Values", "Fitted Values"), col = c("#194f80", "green"), lty = 1, lwd = 2)
plot_HAR_rk_rs
(e) Realized GARCH¶
$ \begin{align} r_t= & \mu_t + \sqrt{h_t}z_t \\ ln(h_t)= & \omega + \beta ln(h_{t-1})+\gamma ln(RV_{t-1})\\ ln(RV_t)= & \xi+ \delta ln(h_t) + u_t \end{align} $
where:
- $z_t$ is a normally distributed iid N(0,1)
- $h_t$ is conditional variance at time ( t )
- $RV_t$ is realized variance at time ( t )
- $ u_t $ is the error term, assumed to be normally distributed.
# Plot the ACF and PACF for the returns
options(repr.plot.width = 10, repr.plot.height = 5)
par(mfrow = c(1, 2))
Acf(amt$ret, main="ACF for Returns")
Pacf(amt$ret, main="PACF for Returns")
auto.arima(amt$ret)
Series: amt$ret
ARIMA(1,0,1) with non-zero mean
Coefficients:
ar1 ma1 mean
0.6467 -0.7596 5e-04
s.e. 0.0950 0.0813 2e-04
sigma^2 = 0.0001936: log likelihood = 4285.44
AIC=-8562.88 AICc=-8562.85 BIC=-8541.62
# Estimate the Realized GARCH model with GARCH order 1,1 and ARMA order 0,0.
real_garchspec= ugarchspec(variance.model = list(model = 'realGARCH', garchOrder = c(1, 1)),
mean.model = list(armaOrder=c(0, 0), include.mean=TRUE), distribution.model = "std")
real_garch_fit= ugarchfit(real_garchspec, amt$ret, solver = 'hybrid', realizedVol = amt$RV)
real_garch_fit
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : realGARCH(1,1)
Mean Model : ARFIMA(0,0,0)
Distribution : std
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.000546 0.000277 1.9698 0.048858
omega 0.774265 0.337744 2.2925 0.021879
alpha1 1.000000 0.154426 6.4756 0.000000
beta1 0.562152 0.056547 9.9413 0.000000
eta11 -0.036534 0.007148 -5.1111 0.000000
eta21 0.026999 0.003208 8.4169 0.000000
delta 0.337887 0.026777 12.6187 0.000000
lambda 0.264883 0.004880 54.2758 0.000000
shape 7.711485 1.321290 5.8363 0.000000
xi -1.653348 0.233836 -7.0705 0.000000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.000546 0.000237 2.2993 0.021488
omega 0.774265 0.487580 1.5880 0.112292
alpha1 1.000000 0.335671 2.9791 0.002891
beta1 0.562152 0.128348 4.3799 0.000012
eta11 -0.036534 0.008461 -4.3177 0.000016
eta21 0.026999 0.004393 6.1465 0.000000
delta 0.337887 0.041407 8.1602 0.000000
lambda 0.264883 0.008095 32.7230 0.000000
shape 7.711485 1.684712 4.5773 0.000005
xi -1.653348 0.359846 -4.5946 0.000004
LogLikelihood : 4356.29
Information Criteria
------------------------------------
Akaike -5.7951
Bayes -5.7596
Shibata -5.7951
Hannan-Quinn -5.7819
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 5.367 0.020516
Lag[2*(p+q)+(p+q)-1][2] 6.294 0.018365
Lag[4*(p+q)+(p+q)-1][5] 10.784 0.005781
d.o.f=0
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 1.396 0.2374
Lag[2*(p+q)+(p+q)-1][5] 1.895 0.6435
Lag[4*(p+q)+(p+q)-1][9] 2.623 0.8193
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 0.2181 0.500 2.000 0.6405
ARCH Lag[5] 0.7950 1.440 1.667 0.7946
ARCH Lag[7] 1.0768 2.315 1.543 0.9007
Nyblom stability test
------------------------------------
Joint Statistic: 3.5441
Individual Statistics:
mu 0.08582
omega 1.16716
alpha1 1.11183
beta1 1.09029
eta11 0.10855
eta21 0.20521
delta 1.54090
lambda 0.07845
shape 0.06165
xi 1.61370
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 2.29 2.54 3.05
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.2121 0.8321
Negative Sign Bias 1.2216 0.2220
Positive Sign Bias 0.5074 0.6119
Joint Effect 2.6308 0.4521
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 24.53 0.1765
2 30 28.00 0.5179
3 40 35.95 0.6099
4 50 50.80 0.4025
Elapsed time : 0.5255949
# Plot acf and pacf of the residuals
res = residuals(real_garch_fit, standardize = TRUE)
sqr_res = residuals(real_garch_fit, standardize = TRUE)^2
options(repr.plot.width = 10, repr.plot.height = 5)
par(mfrow = c(1,2))
Acf(res, lag.max = 20, main="ACF standardized residuals")
Pacf(res, lag.max = 20, main="PACF standardized residuals")
par(mfrow = c(1,2))
Acf(sqr_res, lag.max = 20, main="ACF squared standardized residuals")
Pacf(sqr_res, lag.max = 20, main="PACF squared standardized residuals")
After plotting the residuals we can see that the ACF and PACF still exhibit statistically significant autocorrelation for the standardized residuals, the Weighted Ljung-Box Test on Standardized Residuals in the GARCH output also supports this. Therefore, we decide to reestimate the Realized GARCH model by adding ARMA of order (1,1) as the mean model.
# Re-estimate Realized GARCH with order (1,1) and ARMA order (1,1).
real_garchspec= ugarchspec(variance.model = list(model = 'realGARCH', garchOrder = c(1, 1)),
mean.model = list(armaOrder=c(1, 1), include.mean=TRUE), distribution.model = "std")
real_garch_fit= ugarchfit(real_garchspec, amt$ret, solver = 'hybrid', realizedVol = amt$RV)
real_garch_fit
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : realGARCH(1,1)
Mean Model : ARFIMA(1,0,1)
Distribution : std
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.000614 0.000222 2.7705 0.005597
ar1 0.564222 0.135101 4.1763 0.000030
ma1 -0.650770 0.123916 -5.2517 0.000000
omega 0.733709 0.346462 2.1177 0.034199
alpha1 1.000000 0.156955 6.3712 0.000000
beta1 0.558088 0.056818 9.8224 0.000000
eta11 -0.038538 0.007165 -5.3787 0.000000
eta21 0.026184 0.003154 8.3020 0.000000
delta 0.339895 0.027810 12.2219 0.000000
lambda 0.264784 0.004885 54.2045 0.000000
shape 7.635044 1.283187 5.9501 0.000000
xi -1.632150 0.243184 -6.7116 0.000000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.000614 0.000224 2.7392 0.006159
ar1 0.564222 0.118880 4.7462 0.000002
ma1 -0.650770 0.107331 -6.0632 0.000000
omega 0.733709 0.504072 1.4556 0.145513
alpha1 1.000000 0.341398 2.9291 0.003399
beta1 0.558088 0.128681 4.3370 0.000014
eta11 -0.038538 0.008376 -4.6011 0.000004
eta21 0.026184 0.004192 6.2467 0.000000
delta 0.339895 0.044332 7.6671 0.000000
lambda 0.264784 0.008118 32.6169 0.000000
shape 7.635044 1.655984 4.6106 0.000004
xi -1.632150 0.385393 -4.2350 0.000023
LogLikelihood : 4365.504
Information Criteria
------------------------------------
Akaike -5.8047
Bayes -5.7622
Shibata -5.8048
Hannan-Quinn -5.7888
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.2884 0.5912
Lag[2*(p+q)+(p+q)-1][5] 1.9879 0.9603
Lag[4*(p+q)+(p+q)-1][9] 3.3855 0.8240
d.o.f=2
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 0.736 0.3909
Lag[2*(p+q)+(p+q)-1][5] 1.429 0.7573
Lag[4*(p+q)+(p+q)-1][9] 2.143 0.8872
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 0.07328 0.500 2.000 0.7866
ARCH Lag[5] 0.82538 1.440 1.667 0.7854
ARCH Lag[7] 1.18063 2.315 1.543 0.8825
Nyblom stability test
------------------------------------
Joint Statistic: 3.8048
Individual Statistics:
mu 0.12496
ar1 0.19958
ma1 0.21930
omega 1.18028
alpha1 1.13019
beta1 1.10170
eta11 0.10023
eta21 0.18631
delta 1.60628
lambda 0.07588
shape 0.06057
xi 1.68397
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 2.69 2.96 3.51
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.05354 0.9573
Negative Sign Bias 0.93631 0.3493
Positive Sign Bias 0.29774 0.7659
Joint Effect 1.88631 0.5963
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 20.37 0.3724
2 30 27.12 0.5652
3 40 39.31 0.4561
4 50 54.93 0.2600
Elapsed time : 0.5418198
# Plot acf and pacf of the residuals
res = residuals(real_garch_fit, standardize = TRUE)
sqr_res = residuals(real_garch_fit, standardize = TRUE)^2
options(repr.plot.width = 10, repr.plot.height = 5)
par(mfrow = c(1,2))
Acf(res, lag.max = 20, main="ACF standardized residuals")
Pacf(res, lag.max = 20, main="PACF standardized residuals")
par(mfrow = c(1,2))
Acf(sqr_res, lag.max = 20, main="ACF squared standardized residuals")
Pacf(sqr_res, lag.max = 20, main="PACF squared standardized residuals")
After re-estimating the model and examining the autocorrelations, we can see that the ACF and PACF plots exhibit only minimal autocorrelation in the standarized residuals. Since the autocorrelations at lags 15 and 19 are at the borderline of being statistically significant, we conclude that the model is reasonably well specified, and we would probably not be able to find a parsimonous model that would get rid of these dependencies.
# Plot the true versus fitted values.
amt$real_garch_fit_val = sigma(real_garch_fit)
options(repr.plot.width = 10, repr.plot.height = 5)
plot_real_garch = plot(amt$RV, type = 'l', col = "#194f80", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - realized GARCH")
plot_real_garch = lines(amt$real_garch_fit_val, col = "green", lwd = 1)
plot_real_garch = addLegend("topright", legend.names = c("True Values", "Fitted Values"), col = c("#194f80", "green"), lty = 1, lwd = 2)
plot_real_garch
We can see that the volatility is consistantly overestimated by the realGARCH model.
(f) ARMA-GARCH¶
$ \begin{align} r_t = & \phi_0 + \sum_{i=1}^p \phi_i r_{t-i} + \sum_{j=1}^q \theta_j \epsilon_{t-j} + \epsilon_t, \\ \epsilon_t = & \sigma_t z_t, \quad z_t \sim N(0,1), \\ \sigma_t^2 = & \alpha_0 + \sum_{k=1}^r \alpha_k \epsilon_{t-k}^2 + \sum_{l=1}^s \beta_l \sigma_{t-l}^2, \end{align} $
where:
- $r_t$ are the returns at time ( t )
- $\epsilon_t$ are the residuals at time ( t )
- $\sigma_t^2$ is the conditional variance at time ( t )
- $z_t$ is a iid normally distributed random variable
# Plot the acf and pacf
options(repr.plot.width = 10, repr.plot.height = 5)
par(mfrow = c(1, 2))
Acf(amt$ret, main="ACF for Returns")
Pacf(amt$ret, main="PACF for Returns")
The ACF and PACF plots for the ARMA-GARCH model exhibit autocorrelations at multiple lags. Specifically, we observe significant autocorrelation at lags 1, 3, and 5. Therefore, we will start our analysis with the ARMA(1,1) model specification.
# Let R estimate the model and choose the model specification
auto.arima(amt$ret)
Series: amt$ret
ARIMA(1,0,1) with non-zero mean
Coefficients:
ar1 ma1 mean
0.6467 -0.7596 5e-04
s.e. 0.0950 0.0813 2e-04
sigma^2 = 0.0001936: log likelihood = 4285.44
AIC=-8562.88 AICc=-8562.85 BIC=-8541.62
We let R automatically choose the ARIMA model specification by calling auto.arima. R selected an ARMA(1,1) model, confirming our decision to start with the ARMA(1,1) model specification when estimating the ARMA-GARCH model. The estimated ARMA(1,1) coefficients for the autoregression, moving average, and intercept are statistically significant.
# Estimate the ARMA-GARCH model manually
model_specification = ugarchspec(mean.model = list(armaOrder = c(1, 1)),
variance.model = list(garchOrder = c(1, 1)))
# fit the model
garch_fit = ugarchfit(spec = model_specification, data = amt$ret)
garch_fit
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model : ARFIMA(1,0,1)
Distribution : norm
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.000530 0.000226 2.3437 0.019093
ar1 0.682446 0.107281 6.3613 0.000000
ma1 -0.780022 0.092408 -8.4410 0.000000
omega 0.000005 0.000001 5.3955 0.000000
alpha1 0.054719 0.006214 8.8054 0.000000
beta1 0.920648 0.008867 103.8286 0.000000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.000530 0.000254 2.0854 0.037029
ar1 0.682446 0.139485 4.8926 0.000001
ma1 -0.780022 0.117269 -6.6516 0.000000
omega 0.000005 0.000002 2.3235 0.020152
alpha1 0.054719 0.012435 4.4006 0.000011
beta1 0.920648 0.010465 87.9735 0.000000
LogLikelihood : 4347.196
Information Criteria
------------------------------------
Akaike -5.7883
Bayes -5.7670
Shibata -5.7883
Hannan-Quinn -5.7803
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.2261 0.6344
Lag[2*(p+q)+(p+q)-1][5] 1.4081 0.9991
Lag[4*(p+q)+(p+q)-1][9] 2.9477 0.8984
d.o.f=2
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 6.615 0.01011
Lag[2*(p+q)+(p+q)-1][5] 7.201 0.04634
Lag[4*(p+q)+(p+q)-1][9] 8.147 0.12065
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 0.1144 0.500 2.000 0.7352
ARCH Lag[5] 1.1592 1.440 1.667 0.6863
ARCH Lag[7] 1.5664 2.315 1.543 0.8079
Nyblom stability test
------------------------------------
Joint Statistic: 2.1996
Individual Statistics:
mu 0.07842
ar1 0.77198
ma1 0.73425
omega 0.33583
alpha1 0.25429
beta1 0.25037
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 1.49 1.68 2.12
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.48917 0.6247965
Negative Sign Bias 3.64660 0.0002748 ***
Positive Sign Bias 0.09978 0.9205349
Joint Effect 17.12341 0.0006666 ***
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 80.11 1.783e-09
2 30 105.72 1.165e-10
3 40 108.32 1.925e-08
4 50 125.07 1.409e-08
Elapsed time : 0.09154987
We start estimating the ARMA-GARCH model with an ARMA(1,1) and GARCH(1,1) model specification. All estimates are statistically significant at least at the 5% level. We continue by plotting the ACF and PACF of the standardized residuals to examine how our model fits the data.
res = residuals(garch_fit, standardize = TRUE)
sqr_res = residuals(garch_fit, standardize = TRUE)^2
par(mfrow = c(1,2))
Acf(res, lag.max = 20, main="ACF standardized residuals")
Pacf(res, lag.max = 20, main="PACF standardized residuals")
par(mfrow = c(1,2))
Acf(sqr_res, lag.max = 20, main="ACF squared standardized residuals")
Pacf(sqr_res, lag.max = 20, main="PACF squared standardized residuals")
The ACF and PACF of the standardized residuals show no statistically significant autocorrelation except at lag 19. However, the ACF and PACF of the squared standardized residuals still show statistically significant autocorrelation at lag 1. Therefore, we decided to change the model specification and re-estimate the model.
qqnorm(residuals(garch_fit), cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
qqline(residuals(garch_fit), lwd = 2)
The plot above shows the sample quantiles plotted against the theoretical quantiles. Initially, we assumed that the residuals of our estimated model were normally distributed. However, upon observing the Normal-Q-Q-plot we can see a clear deviation from this assumption, particularly in the tails of our sample distribution. Consequently, we dropped the normality assumption and instead assume that the residuals are t-distributed.
model_specification = ugarchspec(mean.model = list(armaOrder = c(1, 1)),
variance.model = list(garchOrder = c(1, 1)), distribution.model = "std")
# fit the model
garch_t_fit = ugarchfit(spec = model_specification, data = amt$ret)
garch_t_fit
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model : ARFIMA(1,0,1)
Distribution : std
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.000747 0.000214 3.4870 0.000488
ar1 0.665198 0.123016 5.4074 0.000000
ma1 -0.754385 0.108447 -6.9562 0.000000
omega 0.000003 0.000003 1.0348 0.300762
alpha1 0.047463 0.015475 3.0671 0.002162
beta1 0.936158 0.020790 45.0297 0.000000
shape 4.896405 0.605551 8.0859 0.000000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.000747 0.000227 3.28832 0.001008
ar1 0.665198 0.145776 4.56316 0.000005
ma1 -0.754385 0.127814 -5.90219 0.000000
omega 0.000003 0.000013 0.25582 0.798091
alpha1 0.047463 0.046561 1.01938 0.308025
beta1 0.936158 0.073949 12.65956 0.000000
shape 4.896405 0.784868 6.23851 0.000000
LogLikelihood : 4425.839
Information Criteria
------------------------------------
Akaike -5.8918
Bayes -5.8670
Shibata -5.8918
Hannan-Quinn -5.8825
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.03493 0.8517
Lag[2*(p+q)+(p+q)-1][5] 1.36700 0.9994
Lag[4*(p+q)+(p+q)-1][9] 2.94448 0.8989
d.o.f=2
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 8.357 0.003843
Lag[2*(p+q)+(p+q)-1][5] 8.826 0.018310
Lag[4*(p+q)+(p+q)-1][9] 9.772 0.056114
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 0.03122 0.500 2.000 0.8597
ARCH Lag[5] 1.06266 1.440 1.667 0.7144
ARCH Lag[7] 1.48054 2.315 1.543 0.8253
Nyblom stability test
------------------------------------
Joint Statistic: 13.1819
Individual Statistics:
mu 0.1005
ar1 0.5942
ma1 0.6114
omega 1.3286
alpha1 0.2340
beta1 0.1945
shape 0.1283
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 1.69 1.9 2.35
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.63953 0.5225750
Negative Sign Bias 3.84355 0.0001264 ***
Positive Sign Bias 0.06254 0.9501426
Joint Effect 18.16239 0.0004072 ***
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 24.91 0.16364
2 30 46.56 0.02063
3 40 54.40 0.05163
4 50 58.80 0.15936
Elapsed time : 0.1219029
Here, we re-estimate the ARMA-GARCH model with ARMA(1,1) and GARCH(1,1) model specification with assumed Student's t-distributed residuals. Again, all our estimated coefficients, except one, are statistically significant at least at the 5% level. We continue with plotting the residuals to assess how our model fits the data.
res = residuals(garch_t_fit, standardize = TRUE)
sqr_res = residuals(garch_t_fit, standardize = TRUE)^2
par(mfrow = c(1,2))
Acf(res, lag.max = 20, main="ACF standardized residuals")
Pacf(res, lag.max = 20, main="PACF standardized residuals")
par(mfrow = c(1,2))
Acf(sqr_res, lag.max = 20, main="ACF squared standardized residuals")
Pacf(sqr_res, lag.max = 20, main="PACF squared standardized residuals")
After plotting the ACF and PACF of our model's residuals we can still observe persistent autocorrelation of the squared standardized residuals at lag 1.
shape_value = coef(garch_t_fit)["shape"]
print(shape_value)
n_resid = length(residuals(garch_t_fit))
shape 4.896405
The Student's t-distribution has approximately $v$ = 5 degrees of freedom.
qqplot(rt(n_resid,df = shape_value), as.numeric(residuals(garch_t_fit) / sigma(garch_t_fit)), ylab = 'Sample Quantiles',
xlab = 'Theoretical Quantiles', main = 'Student\'s t Q-Q Plot')
qqline(as.numeric(residuals(garch_t_fit)/ sigma(garch_t_fit)))
Here, we plot our sample quantiles against the theoretical quantiles of the Student's t-distribution. We conclude that the assuming Student's t-distributed residuals is a more appropriate assumption than assuming normality. Therefore, we continue estimating our model with assuming Student's t-distributed residuals.
p_max = 3
q_max = 3
ic_min = Inf
best_p = 0
best_q = 0
for (i1 in 1:p_max) {
for (i2 in 1:q_max) {
print(i1, i2)
model_specification = ugarchspec(mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
variance.model = list(garchOrder = c(i1, i2)), distribution.model = "std")
fit = ugarchfit(spec = model_specification, data = amt$ret)
inf_crit = infocriteria(fit)[1] # 1 for aic, 2 for bic, 3 Shibata, 4 Hannan-Quinn
ic_min = ifelse(inf_crit < ic_min, inf_crit, ic_min)
best_p = ifelse(inf_crit == ic_min, i1, best_p)
best_q = ifelse(inf_crit == ic_min, i2, best_q)
}
}
c(best_p, best_q)
[1] 1 [1] 1 [1] 1 [1] 2 [1] 2 [1] 2 [1] 3 [1] 3 [1] 3
- 1
- 3
We iterate over 9 different ARMA(1,1)-GARCH(x,y) model specifications to identify which GARCH specification fits our data best. We assess goodness of fit by using the Akaike Information Criterion (AIC), which evaluates model fit by favoring more parsimonious models. We find that ARMA(0,0)-GARCH(1,3) fits best. Consequently, we proceeded to re-estimate our model using the ARMA(1,1)-GARCH(1,3) specification.
garch_t_1131_spec = ugarchspec(mean.model = list(armaOrder = c(1, 1)),
variance.model = list(model = "sGARCH",garchOrder = c(1, 3)), distribution.model = "std")
# fit the model
garch_t_1131_fit = ugarchfit(spec = garch_t_1131_spec, data = amt$ret)
garch_t_1131_fit
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : sGARCH(1,3)
Mean Model : ARFIMA(1,0,1)
Distribution : std
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.000778 0.000214 3.62969 0.000284
ar1 0.658481 0.129457 5.08647 0.000000
ma1 -0.747423 0.114114 -6.54977 0.000000
omega 0.000007 0.000003 2.29452 0.021761
alpha1 0.105397 0.008582 12.28083 0.000000
beta1 0.249929 0.201909 1.23783 0.215779
beta2 0.148272 0.215624 0.68764 0.491678
beta3 0.461233 0.147228 3.13279 0.001732
shape 5.070949 0.609852 8.31505 0.000000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.000778 0.000220 3.53556 0.000407
ar1 0.658481 0.153519 4.28925 0.000018
ma1 -0.747423 0.135322 -5.52331 0.000000
omega 0.000007 0.000006 1.24485 0.213187
alpha1 0.105397 0.025731 4.09609 0.000042
beta1 0.249929 0.253664 0.98528 0.324489
beta2 0.148272 0.215257 0.68881 0.490941
beta3 0.461233 0.131494 3.50763 0.000452
shape 5.070949 0.692511 7.32255 0.000000
LogLikelihood : 4429.621
Information Criteria
------------------------------------
Akaike -5.8942
Bayes -5.8623
Shibata -5.8942
Hannan-Quinn -5.8823
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.1183 0.7309
Lag[2*(p+q)+(p+q)-1][5] 1.2833 0.9997
Lag[4*(p+q)+(p+q)-1][9] 2.7970 0.9190
d.o.f=2
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 2.794 0.09462
Lag[2*(p+q)+(p+q)-1][11] 4.256 0.69437
Lag[4*(p+q)+(p+q)-1][19] 6.011 0.87291
d.o.f=4
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[5] 1.407 0.500 2.000 0.2355
ARCH Lag[7] 1.676 1.473 1.746 0.5789
ARCH Lag[9] 2.279 2.402 1.619 0.7039
Nyblom stability test
------------------------------------
Joint Statistic: 3.8825
Individual Statistics:
mu 0.1255
ar1 0.5688
ma1 0.5921
omega 0.4559
alpha1 0.2230
beta1 0.1749
beta2 0.1656
beta3 0.1661
shape 0.1178
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 2.1 2.32 2.82
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
t-value prob sig
Sign Bias 0.4510 0.65203
Negative Sign Bias 2.5009 0.01249 **
Positive Sign Bias 0.6581 0.51056
Joint Effect 9.2992 0.02557 **
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 29.95 0.05248
2 30 36.96 0.14725
3 40 51.57 0.08566
4 50 64.07 0.07284
Elapsed time : 0.1368151
Here, we re-estimate the ARMA-GARCH model with ARMA(1,1) and GARCH(1,3) model specification with assumed Student's t-distributed residuals. Many of our estimated coefficients, except three, are statistically significant at least at the 5% level. We continue with plotting the residuals to assess how our model fits the data.
Out of betas, only beta3 is statistically significant, suggesting some persistance of the past conditional volatility. alpha1 is also statistically significant and positive indicating effect of new information.
res = residuals(garch_t_1131_fit, standardize = TRUE)
sqr_res = residuals(garch_t_1131_fit, standardize = TRUE)^2
par(mfrow = c(1,2))
Acf(res, lag.max = 20, main="ACF standardized residuals")
Pacf(res, lag.max = 20, main="PACF standardized residuals")
par(mfrow = c(1,2))
Acf(sqr_res, lag.max = 20, main="ACF squared standardized residuals")
Pacf(sqr_res, lag.max = 20, main="PACF squared standardized residuals")
After plotting the ACF and PACF of our we find no statistically significant autocorrelation at any lag up to lag 20, except at lag 19. We conclude that the ARMA(1,1)-GARCH(1,3) model is reasonably well specified and that increasing the GARCH order helped with modelling the dependencies.
shape_value = coef(garch_t_1131_fit)["shape"]
print(shape_value)
n_resid = length(residuals(garch_t_1131_fit))
shape 5.070949
qqplot(rt(n_resid,df = shape_value), as.numeric(residuals(garch_t_1131_fit) / sigma(garch_t_1131_fit)), ylab = 'Sample Quantiles',
xlab = 'Theoretical Quantiles', main = 'Student\'s t Q-Q Plot')
qqline(as.numeric(residuals(garch_t_1131_fit)/ sigma(garch_t_1131_fit)))
When plotting the sample quantiles against the theoretical quantiles, we can see that the residuals almost perfectly follow the Student's t-distribution. This confirms our choice of using an ARMA(1,1)-GARCH(1,3)-t model.
Box.test(res, type = "Ljung-Box", lag = 25)
Box.test(sqr_res, type = "Ljung-Box", lag = 25)
Box-Ljung test data: res X-squared = 23.394, df = 25, p-value = 0.5545
Box-Ljung test data: sqr_res X-squared = 16.844, df = 25, p-value = 0.8874
We use the Ljung-Box test to test whether our ARMA(1,1)-GARCH(1,3)-t model exhibits any autocorrelation in the residuals. We cannot reject the Null hypothesis of no statistically significant estimate up to the lag 25 and, therefore, conclude that our model caputures the relationships in the data well.
amt$arma_garch_fit_val = sigma(garch_t_1131_fit)
plot_arma_garch = plot(amt$RV, type = 'l', col = "#194f80", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - ARMA-GARCH")
plot_arma_garch = lines(amt$arma_garch_fit_val, col = "green", lwd = 1)
plot_arma_garch = addLegend("topright", legend.names = c("True Values", "Fitted Values"), col = c("#194f80", "green"), lty = 1, lwd = 2)
plot_arma_garch
RV is also seriously overestimated by the ARMA-GARCH, based on the visual inspection it seems that perhaps even the most out of the models we estimated.
(g) Comparison of fits¶
mse_in_samp <- function(true_values, predicted_values) {
mean((true_values - predicted_values)^2, na.rm = TRUE)
}
mae_in_samp <- function(true_values, predicted_values) {
mean(abs(true_values - predicted_values), na.rm = TRUE)
}
fit_cols <- grep("_fit_val", colnames(amt), value = TRUE)
# Calculate MSE and MAE for each model
results <- sapply(fit_cols, function(col) {
fitted_values <- amt[, col]
c(
MSE = mse_in_samp(amt$RV, fitted_values),
MAE = mae_in_samp(amt$RV, fitted_values)
)
})
results_df <- data.frame(t(data.frame(results)))
results_df <- results_df[order(results_df$MSE, results_df$MAE), ] # Sort results by MSE, then by MAE
results_df
| MSE | MAE | |
|---|---|---|
| <dbl> | <dbl> | |
| HAR_rs_rk_fit_val | 1.573535e-05 | 0.002322391 |
| HAR_semi_fit_val | 1.577569e-05 | 0.002335110 |
| HAR_fit_val | 1.606109e-05 | 0.002359915 |
| rv_ar1_fit_val | 1.756837e-05 | 0.002536676 |
| real_garch_fit_val | 2.598096e-05 | 0.003586030 |
| arma_garch_fit_val | 3.091567e-05 | 0.004308790 |
The in-sample MSE and MAE also support the prevous analysis: the GARCH-type models have the highest MAE and MSE.
options(repr.plot.width = 15, repr.plot.height = 10)
par(mfrow = c(3,2))
plot_rv_ar
plot_HAR
plot_HAR_semi
plot_HAR_rk_rs
plot_real_garch
plot_arma_garch
# plot all in one graph
p = plot(amt$RV, type = 'l', col = "#194f80", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - Models")
p = lines(amt$rv_ar1_fit_val, col = "blue", lwd = 1) # AR(1)-RV
p = lines(amt$HAR_fit_val, col = "red", lwd = 1) # HAR
p = lines(amt$HAR_semi_fit_val, col = "green", lwd = 1) # HAR w/ semi-vol
p = lines(amt$HAR_rs_rk_fit_val, col = "purple", lwd = 1) # HAR w/ RK & RS
p = lines(amt$real_garch_fit_val, col = "yellow", lwd = 1) # Realized GARCH
p = lines(amt$arma_garch_fit_val, col = "brown", lwd = 1) # ARMA-GARCH
p = addLegend("topright", legend.names = c("True Values", "AR(1)-RV", "HAR", "HAR w/ Semi-Vol", "HAR w/ RK & RS", "Realized GARCH", "ARMA-GARCH"),
col = c("#194f80", "blue", "red", "green", "purple", "yellow", "brown"), lty = 1, lwd = 2)
p
Here, we plot all fitted values of all estimated models (AR(1)-RV, HAR, HAR with Semi-Volatility, HAR with Realized Kurtosis and Skewness, Realized GARCH, and ARMA-GARCH) against the true values. We find that the Realized GARCH and ARMA-GARCH over-estimate volatility. The other models seem to fit the data well in comparison. Next, we will continue with the forcasting part of our project.
4) Forecasts¶
(a) Running the forecasts¶
# function to perform both expanding and rolling window forecasts for HAR models
forecast_har = function(formula, data, window_length) {
vars = all.vars(formula)
vars = vars[vars != "RV"] # defines the independent varianbles
# here we check where to start - e.g. RV_22 has missing values at the begining of the interval so we have to start later to have 750 obs.
max_non_na_start = max(sapply(data[, vars], function(x) min(which(!is.na(x)))))
start_window = max(window_length, max_non_na_start + 749)
out_of_sample = nrow(data) - start_window
# prepare where to store the results
expanding_forecasts = numeric(out_of_sample)
rolling_forecasts = numeric(out_of_sample)
forecast_dates = index(data)[(start_window + 1):(start_window + out_of_sample)] # for these dates we forecast
for (i in 1:out_of_sample) {
# Expanding window setup
expand_end = start_window + i - 1 # gradually expanding the window from the 750th non-missing obs.
model_expanding = lm(formula, data = data[1:expand_end, ])
new_data_expand = data[expand_end + 1, names(coef(model_expanding))[-1], drop = FALSE]
expanding_forecasts[i] = predict(model_expanding, newdata = new_data_expand)
# Rolling window setup
roll_start = expand_end - window_length + 1 # for roll start is 'window_length' away from the end of the interval
model_rolling = lm(formula, data = data[roll_start:expand_end, ])
new_data_roll = data[expand_end + 1, names(coef(model_rolling))[-1], drop = FALSE]
rolling_forecasts[i] = predict(model_rolling, newdata = new_data_roll)
}
# xts objects for each type of forecast
expanding_xts = xts(expanding_forecasts, order.by = forecast_dates)
rolling_xts = xts(rolling_forecasts, order.by = forecast_dates)
return(list(expanding = expanding_xts, rolling = rolling_xts))
}
# Model specifications
models = list(
RV ~ RV_1,
RV ~ RV_1 + RV_5 + RV_22,
RV ~ RV_p_l + RV_n_l + RV_5 + RV_22,
RV ~ RV_1 + RV_5 + RV_22 + RS_l + RK_l
)
forecast_results = lapply(models, function(model_formula) {
forecast_har(model_formula, amt, 750)
})
forecast_matrix = do.call(cbind, lapply(forecast_results, function(result) { #combine results into a data frame
cbind(result$expanding, result$rolling)}))
# Set column names
colnames(forecast_matrix) = c("AR1_exp", "AR1_roll", "HAR_exp", "HAR_roll", "HAR_semi_exp", "HAR_semi_roll", "HAR_skew_kurt_exp", "HAR_skew_kurt_roll")
head(forecast_matrix) # HAR models start later due to the weekly/montly RV and required 750 obs.
AR1_exp AR1_roll HAR_exp HAR_roll HAR_semi_exp HAR_semi_roll
2013-01-14 0.009430176 0.009430176 NA NA NA NA
2013-01-15 0.009155873 0.009156121 NA NA NA NA
2013-01-16 0.009578290 0.009578171 NA NA NA NA
2013-01-17 0.009332795 0.009335712 NA NA NA NA
2013-01-18 0.008570872 0.008572408 NA NA NA NA
2013-01-22 0.008311133 0.008312711 NA NA NA NA
HAR_skew_kurt_exp HAR_skew_kurt_roll
2013-01-14 NA NA
2013-01-15 NA NA
2013-01-16 NA NA
2013-01-17 NA NA
2013-01-18 NA NA
2013-01-22 NA NA
# function to perform the GARCH forecasting
forecast_garch = function(spec, data, realized_vol, window_length = 750) {
n = nrow(data)
expanding_forecasts = numeric(n - window_length) # objects to store results
rolling_forecasts = numeric(n - window_length)
forecast_dates = index(data)[(window_length + 1):n]
# Expanding window forecasts
for (i in window_length:(n - 1)) {
fit_expanding = ugarchfit(spec = spec, data = data[1:i],
realizedVol = realized_vol[1:i], solver = "hybrid")
forecast_expanding = ugarchforecast(fit_expanding, n.ahead = 1, n.roll = 0)
expanding_forecasts[i - window_length + 1] = forecast_expanding@forecast$sigmaFor[1]}
# Rolling window forecasts
for (i in 1:(n - window_length)) { # starting point from 1 to n - window_length
fit_rolling = ugarchfit(spec = spec, data = data[i:(i + window_length - 1)], # use data from i and the following 750 obs.
realizedVol = realized_vol[i:(i + window_length - 1)], solver = "hybrid")
forecast_rolling = ugarchforecast(fit_rolling, n.ahead = 1, n.roll = 0)
rolling_forecasts[i] = forecast_rolling@forecast$sigmaFor[1]}
expanding_xts = xts(expanding_forecasts, order.by = forecast_dates) # store results
rolling_xts = xts(rolling_forecasts, order.by = forecast_dates)
return(list(expanding = expanding_xts, rolling = rolling_xts))
}
# perform forecasts - specifications defined before
garch_results = list(
realGARCH = forecast_garch(real_garchspec, amt$ret, amt$RV),
sGARCH = forecast_garch(garch_t_1131_spec, amt$ret, amt$RV)
)
# Combine GARCH results analogous to HAR results
garch_forecast_matrix = do.call(cbind, lapply(garch_results, function(result) {
cbind(result$expanding, result$rolling)
}))
# Set column names for GARCH forecasts
colnames(garch_forecast_matrix) = c("realGARCH_exp", "realGARCH_roll", "sGARCH_exp", "sGARCH_roll")
# Combine HAR and GARCH forecasts into one object
combined_forecast_matrix = cbind(forecast_matrix, garch_forecast_matrix)
head(garch_forecast_matrix)
tail(garch_forecast_matrix)
realGARCH_exp realGARCH_roll sGARCH_exp sGARCH_roll 2013-01-11 0.003427967 0.003425051 0.01227074 0.01227074 2013-01-14 0.003467723 0.003448619 0.01204058 0.01201544 2013-01-15 0.003273789 0.003238734 0.01133688 0.01134613 2013-01-16 0.003395814 0.003353558 0.01124976 0.01122078 2013-01-17 0.003381358 0.003365259 0.01152957 0.01148815 2013-01-18 0.002978538 0.002927804 0.01117909 0.01115483
realGARCH_exp realGARCH_roll sGARCH_exp sGARCH_roll 2016-01-14 0.005985543 0.005177259 0.01224169 0.01206505 2016-01-15 0.006718223 0.005706206 0.01441585 0.01388665 2016-01-19 0.009062396 0.007720753 0.01544474 0.01403585 2016-01-20 0.012140318 0.010240255 0.01368778 0.01256856 2016-01-21 0.013886017 0.011727784 0.01362882 0.01338605 2016-01-22 0.013875245 0.011692825 0.01477945 0.01361159
summary(combined_forecast_matrix)
Index AR1_exp AR1_roll HAR_exp
Min. :2013-01-11 Min. :0.007839 Min. :0.007369 Min. :0.006301
1st Qu.:2013-10-10 1st Qu.:0.009331 1st Qu.:0.008939 1st Qu.:0.008528
Median :2014-07-14 Median :0.009948 Median :0.009509 Median :0.009418
Mean :2014-07-14 Mean :0.010270 Mean :0.009833 Mean :0.009629
3rd Qu.:2015-04-14 3rd Qu.:0.010742 3rd Qu.:0.010264 3rd Qu.:0.010260
Max. :2016-01-22 Max. :0.031653 Max. :0.031704 Max. :0.022883
NA's :1 NA's :1 NA's :22
HAR_roll HAR_semi_exp HAR_semi_roll HAR_skew_kurt_exp
Min. :0.006129 Min. :0.006245 Min. :0.005892 Min. :0.005322
1st Qu.:0.008492 1st Qu.:0.008502 1st Qu.:0.008452 1st Qu.:0.008513
Median :0.009292 Median :0.009361 Median :0.009222 Median :0.009459
Mean :0.009496 Mean :0.009638 Mean :0.009513 Mean :0.009665
3rd Qu.:0.010065 3rd Qu.:0.010302 3rd Qu.:0.010139 3rd Qu.:0.010374
Max. :0.022672 Max. :0.019494 Max. :0.019161 Max. :0.019949
NA's :22 NA's :22 NA's :22 NA's :22
HAR_skew_kurt_roll realGARCH_exp realGARCH_roll sGARCH_exp
Min. :0.004940 Min. :0.001933 Min. :0.001992 Min. :0.008366
1st Qu.:0.008488 1st Qu.:0.003712 1st Qu.:0.003133 1st Qu.:0.010776
Median :0.009298 Median :0.004511 Median :0.003887 Median :0.012155
Mean :0.009515 Mean :0.004730 Mean :0.005620 Mean :0.012811
3rd Qu.:0.010111 3rd Qu.:0.005295 3rd Qu.:0.004846 3rd Qu.:0.013966
Max. :0.019653 Max. :0.013886 Max. :0.590552 Max. :0.031098
NA's :22
sGARCH_roll
Min. :0.007668
1st Qu.:0.010439
Median :0.011796
Mean :0.012451
3rd Qu.:0.013698
Max. :0.032497
We can see that the forecasts range from about 0.002 to 0.03 for all models and forecasting schemes. The only exception is the realGARCH model forecasted with the rolling window. There appear to be some outliers reaching almost 0.6.
options(repr.plot.width = 5, repr.plot.height = 5)
boxplot(combined_forecast_matrix$realGARCH_exp)
boxplot(combined_forecast_matrix$realGARCH_roll)
Indeed, the boxplot supports that that the presence of outliers.
(b) Forecasts evaluation - graphs¶
# combine data for plotting
data_all = cbind(combined_forecast_matrix, amt$RV)
data_all_nona = na.omit(data_all)
options(repr.plot.width = 10, repr.plot.height = 10)
p = plot(data_all_nona$RV, type = 'l', col = "#194f80", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - Models")
series_list = list(
AR1_exp = "blue",AR1_roll = "darkblue",HAR_exp = "red",
HAR_roll = "darkred",HAR_semi_exp = "green",HAR_semi_roll = "darkgreen",
HAR_skew_kurt_exp = "purple",HAR_skew_kurt_roll = "darkmagenta",
realGARCH_exp = "yellow",realGARCH_roll = "gold",sGARCH_exp = "orange",sGARCH_roll = "darkorange")
# for loop through the list
for (series_name in names(series_list)) {
p = lines(data_all_nona[,series_name], col = series_list[[series_name]], lwd = 1)
}
p = addLegend("topright", legend.names = c("True Values", names(series_list)),col = c("#194f80", unlist(series_list)),lty = 1,lwd = 2)
p
Clearly, the realGARCH forecast produces a few extreme values using the rolling window. Let's also plot the series without it for readability.
options(repr.plot.width = 10, repr.plot.height = 10)
# same code just commenting out realGARCH_roll
p = plot(data_all_nona$RV, type = 'l', col = "black", lwd = 2, xlab = "Time", ylab = "RV", main = "True vs Fitted Values - Models")
series_list = list(
AR1_exp = "blue",
AR1_roll = "darkblue",
HAR_exp = "red",
HAR_roll = "darkred",
HAR_semi_exp = "green",
HAR_semi_roll = "darkgreen",
HAR_skew_kurt_exp = "purple",
HAR_skew_kurt_roll = "darkmagenta",
realGARCH_exp = "gold",
#realGARCH_roll = "yellow",
sGARCH_exp = "orange",
sGARCH_roll = "red"
)
for (series_name in names(series_list)) {
p = lines(data_all_nona[,series_name], col = series_list[[series_name]], lwd = 1)
}
p = addLegend("topright",legend.names = c("True Values", names(series_list)),col = c("black", unlist(series_list)), lty = 1, lwd = 2)
p
While the realGARCH using the expanding window systematically underforecasts the volatility, the standard ARMA-GARCH overforecasts it, regardless of the forecasting schemes.
# calculate errors
error_list = list()
forecast_names = colnames(combined_forecast_matrix)
for (name in forecast_names) {
error_list[[name]] = amt$RV - combined_forecast_matrix[,name]
}
# Plot the errors
options(repr.plot.width = 15, repr.plot.height = 15)
par(mfrow=c(4, 3))
for (i in 1:length(forecast_names)) {
plot(index(combined_forecast_matrix), error_list[[forecast_names[i]]], type='l', col="#194f80", lwd=1,
xlab="Time", ylab="Error", main=paste("Error -", forecast_names[i]))
}
unique(sub("_(exp|roll)$", "", forecast_names)) # how to get the model names
- 'AR1'
- 'HAR'
- 'HAR_semi'
- 'HAR_skew_kurt'
- 'realGARCH'
- 'sGARCH'
# plot roll vs expanding window for each model
par(mfrow=c(3, 2))
for (model in unique(sub("_(exp|roll)$", "", forecast_names))) {
exp_name = paste0(model, "_exp")
roll_name = paste0(model, "_roll")
plot(index(combined_forecast_matrix), error_list[[exp_name]], type='l', col="#194f80", lwd=1,
xlab="Time", ylab="Error", main=paste("Error -", model, "Exp vs Roll"))
lines(index(combined_forecast_matrix), error_list[[roll_name]], col="green", lwd=1)
abline(h = 0, col = "red", lwd = 1)
legend("topright", legend=c("Expanding", "Rolling"), col=c("#194f80", "green"), lty=1, lwd=2)
}
The errors are quite similar across models and when comparing the two forecasting schemes. The most obvious difference is the realGARCH model using the rolling window, which leads to a few outliers in the forecast. Since forecasts with large errors can lead to serious consequences in the real world, we continue the analysis with those forecasts as they are, as well as winsorized values, to see the performance of the model without being considerably affected by the outliers.
# Winsorizing realGARCH_exp at the 99.5th percentile
combined_forecast_matrix$realGARCH_wins_exp = Winsorize(combined_forecast_matrix$realGARCH_exp, probs = c(0.005, 0.995))
# Winsorizing realGARCH_roll at the 99th percentile
combined_forecast_matrix$realGARCH_wins_roll = Winsorize(combined_forecast_matrix$realGARCH_roll, probs = c(0.005, 0.995))
summary(combined_forecast_matrix)
Index AR1_exp AR1_roll HAR_exp
Min. :2013-01-11 Min. :0.007839 Min. :0.007369 Min. :0.006301
1st Qu.:2013-10-10 1st Qu.:0.009331 1st Qu.:0.008939 1st Qu.:0.008528
Median :2014-07-14 Median :0.009948 Median :0.009509 Median :0.009418
Mean :2014-07-14 Mean :0.010270 Mean :0.009833 Mean :0.009629
3rd Qu.:2015-04-14 3rd Qu.:0.010742 3rd Qu.:0.010264 3rd Qu.:0.010260
Max. :2016-01-22 Max. :0.031653 Max. :0.031704 Max. :0.022883
NA's :1 NA's :1 NA's :22
HAR_roll HAR_semi_exp HAR_semi_roll HAR_skew_kurt_exp
Min. :0.006129 Min. :0.006245 Min. :0.005892 Min. :0.005322
1st Qu.:0.008492 1st Qu.:0.008502 1st Qu.:0.008452 1st Qu.:0.008513
Median :0.009292 Median :0.009361 Median :0.009222 Median :0.009459
Mean :0.009496 Mean :0.009638 Mean :0.009513 Mean :0.009665
3rd Qu.:0.010065 3rd Qu.:0.010302 3rd Qu.:0.010139 3rd Qu.:0.010374
Max. :0.022672 Max. :0.019494 Max. :0.019161 Max. :0.019949
NA's :22 NA's :22 NA's :22 NA's :22
HAR_skew_kurt_roll realGARCH_exp realGARCH_roll sGARCH_exp
Min. :0.004940 Min. :0.001933 Min. :0.001992 Min. :0.008366
1st Qu.:0.008488 1st Qu.:0.003712 1st Qu.:0.003133 1st Qu.:0.010776
Median :0.009298 Median :0.004511 Median :0.003887 Median :0.012155
Mean :0.009515 Mean :0.004730 Mean :0.005620 Mean :0.012811
3rd Qu.:0.010111 3rd Qu.:0.005295 3rd Qu.:0.004846 3rd Qu.:0.013966
Max. :0.019653 Max. :0.013886 Max. :0.590552 Max. :0.031098
NA's :22
sGARCH_roll realGARCH_wins_exp realGARCH_wins_roll
Min. :0.007668 Min. :0.002287 Min. :0.002216
1st Qu.:0.010439 1st Qu.:0.003712 1st Qu.:0.003133
Median :0.011796 Median :0.004511 Median :0.003887
Mean :0.012451 Mean :0.004725 Mean :0.004238
3rd Qu.:0.013698 3rd Qu.:0.005295 3rd Qu.:0.004846
Max. :0.032497 Max. :0.012124 Max. :0.011772
# calculate errors
error_list = list()
forecast_names = colnames(combined_forecast_matrix)
for (name in forecast_names) {
error_list[[name]] = amt$RV - combined_forecast_matrix[,name]
colnames(error_list[[name]]) = name
}
# Plot the errors
options(repr.plot.width = 15, repr.plot.height = 15)
par(mfrow=c(5, 3))
for (i in 1:length(forecast_names)) {
plot(index(combined_forecast_matrix), error_list[[forecast_names[i]]], type='l', col="#194f80", lwd=1,
xlab="Time", ylab="Error", main=paste("Error -", forecast_names[i]))
}
model_types = unique(sub("_(exp|roll)$", "", forecast_names))
options(repr.plot.width = 15, repr.plot.height = 15)
par(mfrow=c(3, 2))
for (model in model_types) {
exp_name = paste0(model, "_exp")
roll_name = paste0(model, "_roll")
if (exp_name %in% forecast_names && roll_name %in% forecast_names) {
plot(index(combined_forecast_matrix), error_list[[exp_name]], type='l', col="#194f80", lwd=1,
xlab="Time", ylab="Error", main=paste("Error -", model, "Exp vs Roll"))
lines(index(combined_forecast_matrix), error_list[[roll_name]], col="green", lwd=1)
abline(h = 0, col = "red", lwd = 1)
legend("topright", legend=c("Expanding", "Rolling"), col=c("#194f80", "green"), lty=1, lwd=2)
}
}
(c) Forecasts evaluation - MAE and MSE¶
# Function to calculate MSE
mse = function(errors) {
round(mean(na.omit(errors)^2), 7)
}
# Function to calculate MAE
mae = function(errors) {
round(mean(abs(na.omit(errors))), 7)
}
mse_results = sapply(error_list, mse)
mae_results = sapply(error_list, mae)
print(sort(mae_results, decreasing = FALSE))
HAR_skew_kurt_roll HAR_semi_roll HAR_skew_kurt_exp HAR_roll
0.0021174 0.0021212 0.0021336 0.0021347
HAR_semi_exp HAR_exp AR1_roll AR1_exp
0.0021382 0.0021482 0.0023006 0.0024359
sGARCH_roll sGARCH_exp realGARCH_exp realGARCH_wins_exp
0.0041189 0.0043761 0.0046623 0.0046650
realGARCH_wins_roll realGARCH_roll
0.0051633 0.0065451
print(sort(mse_results, decreasing = FALSE))
HAR_semi_exp HAR_semi_roll HAR_skew_kurt_exp HAR_exp
0.0000125 0.0000126 0.0000127 0.0000128
HAR_skew_kurt_roll HAR_roll AR1_roll AR1_exp
0.0000128 0.0000129 0.0000135 0.0000138
sGARCH_roll sGARCH_exp realGARCH_exp realGARCH_wins_exp
0.0000281 0.0000300 0.0000336 0.0000337
realGARCH_wins_roll realGARCH_roll
0.0000390 0.0007684
The GARCH-type models are outperformed by the heterogeneous autoregressive models in terms of the mean absolute error and the mean squared error. The MAE and MSE is similar for all HAR model, followed by the simple AR(1)-RV model, and then followed by the ARMA-GARCH, which has the MAE and MSE approximately twice as large as as the AR1 model, regardless of the forecasting schemes. Overall, the difference between the rolling and expanding scheme MAE and MSE is relatively minor for all models with the exception of the realGARCH, which is mostly affected by the extreme values discussed before.
(d) Diebold-Mariano test¶
Since we have 12 different models (+ realGARCH with winsorized forecasts), the most systematic way to compare them using the Diebold-Mariano test would be the number of times the certain models outperforms another.
mse_daily = function(errors) {
errors^2 # Squaring each error
}
daily_mse_results = lapply(error_list, mse_daily)
daily_mse_matrix = do.call(cbind, daily_mse_results)
head(daily_mse_matrix)
AR1_exp AR1_roll HAR_exp HAR_roll HAR_semi_exp
2013-01-11 NA NA NA NA NA
2013-01-14 1.119151e-05 1.119151e-05 NA NA NA
2013-01-15 4.484798e-06 4.485850e-06 NA NA NA
2013-01-16 9.420028e-06 9.419292e-06 NA NA NA
2013-01-17 2.012151e-05 2.014768e-05 NA NA NA
2013-01-18 1.817745e-05 1.819055e-05 NA NA NA
HAR_semi_roll HAR_skew_kurt_exp HAR_skew_kurt_roll realGARCH_exp
2013-01-11 NA NA NA 1.055429e-05
2013-01-14 NA NA NA 6.849116e-06
2013-01-15 NA NA NA 1.417033e-05
2013-01-16 NA NA NA 9.692452e-06
2013-01-17 NA NA NA 2.148385e-06
2013-01-18 NA NA NA 1.765794e-06
realGARCH_roll sGARCH_exp sGARCH_roll realGARCH_wins_exp
2013-01-11 1.057324e-05 3.129320e-05 3.129320e-05 1.055429e-05
2013-01-14 6.949473e-06 3.547124e-05 3.517241e-05 6.849116e-06
2013-01-15 1.443548e-05 1.847921e-05 1.855876e-05 1.417033e-05
2013-01-16 9.957342e-06 2.247404e-05 2.220008e-05 9.692452e-06
2013-01-17 2.195837e-06 4.465548e-05 4.410356e-05 2.148385e-06
2013-01-18 1.903202e-06 4.722055e-05 4.688772e-05 1.765794e-06
realGARCH_wins_roll
2013-01-11 1.057324e-05
2013-01-14 6.949473e-06
2013-01-15 1.443548e-05
2013-01-16 9.957342e-06
2013-01-17 2.195837e-06
2013-01-18 1.903202e-06
mae_daily = function(errors) {
abs(errors)
}
daily_mae_results = lapply(error_list, mae_daily)
daily_mae_matrix = do.call(cbind, daily_mse_results)
perform_dm_tests = function(data_matrix, cols_to_remove = NULL) {
data_matrix = na.omit(data_matrix[, -cols_to_remove])
n_models = ncol(data_matrix)
results = matrix(NA, nrow = n_models, ncol = n_models, dimnames = list(colnames(data_matrix), colnames(data_matrix)))
outperformance_counts = integer(n_models)
outperformance_list = vector("list", n_models)
names(outperformance_list) = colnames(data_matrix)
for (i in 1:n_models) {
for (j in 1:n_models) {
if (i != j) {
dm_result = dm.test(data_matrix[, i], data_matrix[, j], alternative = "less")
results[j, i] = dm_result$p.value
if (dm_result$p.value < 0.05) {
outperformance_counts[i] = outperformance_counts[i] + 1
outperformance_list[[i]] = c(outperformance_list[[i]], colnames(data_matrix)[j])}}}}
names(outperformance_counts) = colnames(data_matrix)
sorted_counts = sort(outperformance_counts, decreasing = TRUE)
print("Sorted counts of model outperformance:")
print(sorted_counts)
# Print models each model outperforms
print("Models each model outperforms:")
for (i in names(outperformance_list)) {
if (length(outperformance_list[[i]]) > 0) {
cat(i, "outperforms:", paste(outperformance_list[[i]], collapse=", "), "\n")
}
}
}
perform_dm_tests(daily_mse_matrix, c((ncol(daily_mse_matrix)-1):ncol(daily_mse_matrix)))
[1] "Sorted counts of model outperformance:"
AR1_exp HAR_exp HAR_semi_exp HAR_skew_kurt_exp
2 2 2 2
AR1_roll HAR_roll HAR_semi_roll HAR_skew_kurt_roll
1 1 1 1
realGARCH_exp realGARCH_roll sGARCH_exp sGARCH_roll
0 0 0 0
[1] "Models each model outperforms:"
AR1_exp outperforms: AR1_roll, realGARCH_exp
AR1_roll outperforms: realGARCH_exp
HAR_exp outperforms: HAR_roll, realGARCH_exp
HAR_roll outperforms: realGARCH_exp
HAR_semi_exp outperforms: HAR_semi_roll, realGARCH_exp
HAR_semi_roll outperforms: realGARCH_exp
HAR_skew_kurt_exp outperforms: HAR_skew_kurt_roll, realGARCH_exp
HAR_skew_kurt_roll outperforms: realGARCH_exp
perform_dm_tests(daily_mse_matrix, c(9, 10))
[1] "Sorted counts of model outperformance:"
AR1_exp HAR_exp HAR_semi_exp HAR_skew_kurt_exp
3 3 3 3
AR1_roll HAR_roll HAR_semi_roll HAR_skew_kurt_roll
2 2 2 2
sGARCH_exp sGARCH_roll realGARCH_wins_exp realGARCH_wins_roll
1 1 1 0
[1] "Models each model outperforms:"
AR1_exp outperforms: AR1_roll, realGARCH_wins_exp, realGARCH_wins_roll
AR1_roll outperforms: realGARCH_wins_exp, realGARCH_wins_roll
HAR_exp outperforms: HAR_roll, realGARCH_wins_exp, realGARCH_wins_roll
HAR_roll outperforms: realGARCH_wins_exp, realGARCH_wins_roll
HAR_semi_exp outperforms: HAR_semi_roll, realGARCH_wins_exp, realGARCH_wins_roll
HAR_semi_roll outperforms: realGARCH_wins_exp, realGARCH_wins_roll
HAR_skew_kurt_exp outperforms: HAR_skew_kurt_roll, realGARCH_wins_exp, realGARCH_wins_roll
HAR_skew_kurt_roll outperforms: realGARCH_wins_exp, realGARCH_wins_roll
sGARCH_exp outperforms: realGARCH_wins_roll
sGARCH_roll outperforms: realGARCH_wins_roll
realGARCH_wins_exp outperforms: realGARCH_wins_roll
perform_dm_tests(daily_mae_matrix, c((ncol(daily_mae_matrix)-1):ncol(daily_mae_matrix)))
[1] "Sorted counts of model outperformance:"
AR1_exp HAR_exp HAR_semi_exp HAR_skew_kurt_exp
2 2 2 2
AR1_roll HAR_roll HAR_semi_roll HAR_skew_kurt_roll
1 1 1 1
realGARCH_exp realGARCH_roll sGARCH_exp sGARCH_roll
0 0 0 0
[1] "Models each model outperforms:"
AR1_exp outperforms: AR1_roll, realGARCH_exp
AR1_roll outperforms: realGARCH_exp
HAR_exp outperforms: HAR_roll, realGARCH_exp
HAR_roll outperforms: realGARCH_exp
HAR_semi_exp outperforms: HAR_semi_roll, realGARCH_exp
HAR_semi_roll outperforms: realGARCH_exp
HAR_skew_kurt_exp outperforms: HAR_skew_kurt_roll, realGARCH_exp
HAR_skew_kurt_roll outperforms: realGARCH_exp
perform_dm_tests(daily_mae_matrix, c(9, 10))
[1] "Sorted counts of model outperformance:"
AR1_exp HAR_exp HAR_semi_exp HAR_skew_kurt_exp
3 3 3 3
AR1_roll HAR_roll HAR_semi_roll HAR_skew_kurt_roll
2 2 2 2
sGARCH_exp sGARCH_roll realGARCH_wins_exp realGARCH_wins_roll
1 1 1 0
[1] "Models each model outperforms:"
AR1_exp outperforms: AR1_roll, realGARCH_wins_exp, realGARCH_wins_roll
AR1_roll outperforms: realGARCH_wins_exp, realGARCH_wins_roll
HAR_exp outperforms: HAR_roll, realGARCH_wins_exp, realGARCH_wins_roll
HAR_roll outperforms: realGARCH_wins_exp, realGARCH_wins_roll
HAR_semi_exp outperforms: HAR_semi_roll, realGARCH_wins_exp, realGARCH_wins_roll
HAR_semi_roll outperforms: realGARCH_wins_exp, realGARCH_wins_roll
HAR_skew_kurt_exp outperforms: HAR_skew_kurt_roll, realGARCH_wins_exp, realGARCH_wins_roll
HAR_skew_kurt_roll outperforms: realGARCH_wins_exp, realGARCH_wins_roll
sGARCH_exp outperforms: realGARCH_wins_roll
sGARCH_roll outperforms: realGARCH_wins_roll
realGARCH_wins_exp outperforms: realGARCH_wins_roll
There does not seem to be clear-cut winner in our data in terms of forecasting performance. The general conclusion is that the (H)AR models outperform the GARCH-type models. According to the dm.test it also seems that the expanding window forecast perform better to their rolling window counterparts.
(e) Minzer-Zarnowitz regression¶
Where:
- $ y_t $ : Actual value at time $ t $
- $ \hat{y}_t $ : Forecasted value at time $ t $
- $ \beta_0 $ : Intercept term of the regression
- $ \beta_1 $ : Slope coefficient, which measures the bias in the forecast
- $ \epsilon_t $ : Error term, capturing the unpredicted part of $ y_t $
In the Minzer-Zarnowitz regression we regress the actual values ($y_t$) on the forecasted values ($\hat{y}_t$). We conduct the following hypothesis tests:
Intercept Term ($\beta_0$): We test whether the intercept is statistically different from zero. A non-zero intercept would indicate a systematic bias in the forecast, suggesting that the forecasts are either consistently overestimating or underestimating the actual values.
Slope Coefficient ($\beta_1$): We test whether the slope is statistically different from one. The slope being significantly different from one indicates that the forecasts are not perfectly proportional to the actual values
Joint Hypothesis Test via Wald Test: We use the Wald test to jointly assess the null hypotheses that $\beta_0 = 0$ and $\beta_1 = 1$. This test provides a robust measure of overall model fit and the adequacy of forecasts. If the null hypothesis is rejected, it suggests that the forecasts do not adequately follow the actual values
# empty list to store results
results = list()
# Loop through each predictor column
for (column in names(data_all[, -ncol(data_all)])) {
# estimate model
model = lm(paste("RV ~", column), data = data_all)
# Extract coefficients
coef_info = summary(model)$coefficients
vcov_matrix = vcovHAC(model) # Get the variance-covariance matrix of the coefficients
# Calculate Wald test statistic
wald_statistic = (coef_info[, "Estimate"] - c(0, 1)) %*% solve(vcov_matrix) %*% (coef_info[, "Estimate"] - c(0, 1))
wald_p_value = pchisq(wald_statistic, df = 2, lower.tail = FALSE) # Chi-squared test with 2 degrees of freedom
# Calculate p-values for each coefficient individually
intercept_p_value = 2 * (1 - pnorm(abs(coef_info["(Intercept)", "Estimate"] / coef_info["(Intercept)", "Std. Error"])))
slope_p_value = 2 * (1 - pnorm(abs((coef_info[column, "Estimate"] - 1) / coef_info[column, "Std. Error"])))
# Store results in the list
results[[column]] = c(
round(coef_info[column, "Estimate"], 2),
round(coef_info[column, "Std. Error"], 2),
round((coef_info[column, "Estimate"] - 1) / coef_info[column, "Std. Error"], 2),
round(slope_p_value, 2),
round(intercept_p_value, 2),
round(wald_p_value, 2),
ifelse(slope_p_value < 0.05, "Reject", "Do not reject"), # Reject/Not Reject for slope = 1
ifelse(intercept_p_value < 0.05, "Reject", "Do not reject"), # Reject/Not Reject for intercept = 0
ifelse(wald_p_value < 0.05, "Reject", "Do not reject") # Reject/Not Reject for Wald test
)}
# Convert results to a data frame
results_df = do.call(rbind, results)
colnames(results_df) = c(
"Slope Coeff.",
"SE",
"Predictor Test Stat.",
"Slope p-value",
"Intercept p-value",
"Wald p-value",
"beta==1",
"(Intercept)==0",
"beta==1∧(Intercept)==0"
)
results_df
| Slope Coeff. | SE | Predictor Test Stat. | Slope p-value | Intercept p-value | Wald p-value | beta==1 | (Intercept)==0 | beta==1∧(Intercept)==0 | |
|---|---|---|---|---|---|---|---|---|---|
| AR1_exp | 0.83 | 0.08 | -2.26 | 0.02 | 0.28 | 0 | Reject | Do not reject | Reject |
| AR1_roll | 0.82 | 0.08 | -2.08 | 0.04 | 0.13 | 0 | Reject | Do not reject | Reject |
| HAR_exp | 0.9 | 0.08 | -1.28 | 0.2 | 0.31 | 0.35 | Do not reject | Do not reject | Do not reject |
| HAR_roll | 0.92 | 0.08 | -0.95 | 0.34 | 0.39 | 0.85 | Do not reject | Do not reject | Do not reject |
| HAR_semi_exp | 0.91 | 0.07 | -1.25 | 0.21 | 0.34 | 0.32 | Do not reject | Do not reject | Do not reject |
| HAR_semi_roll | 0.93 | 0.08 | -0.92 | 0.36 | 0.42 | 0.8 | Do not reject | Do not reject | Do not reject |
| HAR_skew_kurt_exp | 0.89 | 0.07 | -1.52 | 0.13 | 0.23 | 0.23 | Do not reject | Do not reject | Do not reject |
| HAR_skew_kurt_roll | 0.91 | 0.08 | -1.09 | 0.28 | 0.33 | 0.76 | Do not reject | Do not reject | Do not reject |
| realGARCH_exp | 1.06 | 0.08 | 0.7 | 0.48 | 0 | 0 | Do not reject | Reject | Reject |
| realGARCH_roll | 0 | 0.01 | -192.64 | 0 | 0 | 0 | Reject | Reject | Reject |
| sGARCH_exp | 0.32 | 0.05 | -14.76 | 0 | 0 | 0 | Reject | Reject | Reject |
| sGARCH_roll | 0.3 | 0.05 | -15.4 | 0 | 0 | 0 | Reject | Reject | Reject |
We do not reject the joint null hypothesis for all heterogenous autoregressive type models, while it is rejected for the AR and GARCH-type models.
5) Summary¶
This project analyzes and forecasts the volatility of the stock price of American Tower Corporation (AMT) using various models and forecasting schemes. We estimate the AR(1)-RV, HAR, HAR-RS, HAR-Rskew-RKurt, Realized GARCH, ARMA-GARCH (specifically ARMA(1,1)-GARCH(1,3)-t) models Then we perform an out-of-sample forecast using rolling and expanding window forecasting schemes.
Observing the in-sample fit, the HAR models perform the best in terms of the MAE and MSE metrics, while the ARMA-GARCH model performs the worst. This can also be seen from the in-sample plots of the fitted values against the true values as the ARMA-GARCH appears to be significantly overestimating the volatility. The HAR models all appear to be performing very similarly; however, the HAR-Rskew-Rkurt does appear to have the smallest in-sample MAE and MSE by a small margin.
The forecasts using the real-GARCH model produce some extreme values for both forecasting schemes. Visual inspection indicates that the ARMA-GARCH model overestimates the volatility, while the realGARCH forecast underforecasts it. The rolling and expanding schemes generally produce similar results with only minor differences. Computing the MAE and MSE of the forecasts, we see that all the HAR models exhibit very similar results and the smallest error metrics, while the real-GARCH has the higher MAE and MSE.
From the Diebold-Mariano test, we can see that forecasts made using the expanding window usually outperform forecasts using the rolling window for the same models, perhaps due to the inclusion of additional observations, and hence the ability to exploit some long-term information which the rolling window could miss. The test also support that the HAR-type models and the AR(1) model perform better than the GARCH-type models.
We use the Minzer-Zarnowitz regression to evaluate our forecast and test formally whether they systematically overestimate or underestimate the volatility and also to uncover potential bias. We do not reject the null hypothesis of $\alpha = 0$ and $\beta = 1$ for all the HAR models.
In conclusion, the HAR models are preferred as they perform better than the GARCH models in every metric. Namely, we would prefer either the HAR-RS or HAR-Rskew-RKurt due to their ability to account for the additional varialbes.